﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// This structure represents the type over which the <see cref="NormalGamma"/> distribution
    /// is defined.
    /// </summary>
    public struct MeanPrecisionPair
    {
        private double mMean;
        private double mPrecision;

        /// <summary>
        /// Constructs a new mean precision pair.
        /// </summary>
        /// <param name="m">The mean of the pair.</param>
        /// <param name="p">The precision of the pair.</param>
        public MeanPrecisionPair(double m, double p)
        {
            mMean = m;
            mPrecision = p;
        }

        /// <summary>
        /// Gets/sets the mean of the pair.
        /// </summary>
        public double Mean
        {
            get { return mMean; }
            set { mMean = value; }
        }

        /// <summary>
        /// Gets/sets the precision of the pair.
        /// </summary>
        public double Precision
        {
            get { return mPrecision; }
            set { mPrecision = value; }
        }
    }

    /// <summary>
    /// The NormalGamma distribution is the conjugate prior distribution for the <see cref="Normal"/>
    /// distribution. It specifies a prior over the mean and precision of the <see cref="Normal"/> distribution.
    /// 
    /// It is parameterized by four numbers: the mean location, the mean scale, the precision shape and the
    /// precision inverse scale. The following degenerate cases are special: when the precision is known,
    /// the precision shape will encode the value of the precision while the precision inverse scale is positive
    /// infinity. When the mean is known, the mean location will encode the value of the mean while the scale
    /// will be 0. A completely degenerate NormalGamma distribution with known mean and precision is possible as well.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class NormalGamma
    {
        // The location of the mean.
        private readonly double mMeanLocation;

        // The scale of the mean.
        private readonly double mMeanScale;

        // The shape of the precision.
        private readonly double mPrecisionShape;

        // The inverse scale of the precision.
        private readonly double mPrecisionInvScale;

        /// <summary>
        /// Constructs a NormalGamma distribution.
        /// </summary>
        /// <param name="meanLocation">The location of the mean.</param>
        /// <param name="meanScale">The scale of the mean.</param>
        /// <param name="precShape">The shape of the precision.</param>
        /// <param name="precInvScale">The inverse scale of the precision.</param>
        public NormalGamma(double meanLocation, double meanScale, double precShape, double precInvScale)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(meanScale, precShape, precInvScale);
            }

            mMeanLocation = meanLocation;
            mMeanScale = meanScale;
            mPrecisionShape = precShape;
            mPrecisionInvScale = precInvScale;
            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        public override string ToString()
        {
            return "NormalGamma(Mean Location = " + mMeanLocation + ", Mean Scale = " + mMeanScale +
                        ", Precision Shape = " + mPrecisionShape + ", Precision Inverse Scale = " + mPrecisionInvScale + ")";
        }

        /// <summary>
        /// Gets the location of the mean.
        /// </summary>
        public double MeanLocation
        {
            get { return mMeanLocation; }
        }

        /// <summary>
        /// Gets the scale of the mean.
        /// </summary>
        public double MeanScale
        {
            get { return mMeanScale; }
        }

        /// <summary>
        /// Gets the shape of the precision.
        /// </summary>
        public double PrecisionShape
        {
            get { return mPrecisionShape; }
        }

        /// <summary>
        /// Gets the inverse scale of the precision.
        /// </summary>
        public double PrecisionInverseScale
        {
            get { return mPrecisionInvScale; }
        }

        /// <summary>
        /// Gets the mean of the distribution.
        /// </summary>
        /// <value>The mean of the distribution.</value>
        public MeanPrecisionPair Mean
        {
            get
            {
                if (Double.IsPositiveInfinity(mPrecisionInvScale))
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape);
                }
                else
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape / mPrecisionInvScale);
                }
            }
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        /// <summary>
        /// The mode of the distribution.
        /// </summary>
        /// <value></value>
        public MeanPrecisionPair Mode
        {
            get
            {
                if (Double.IsPositiveInfinity(mPrecisionInvScale))
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape);
                }
                else
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape / mPrecisionInvScale);
                }
            }
        }

        /// <summary>
        /// The median of the distribution.
        /// </summary>
        /// <value></value>
        public MeanPrecisionPair Median
        {
            get
            {
                if (Double.IsPositiveInfinity(mPrecisionInvScale))
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape);
                }
                else
                {
                    return new MeanPrecisionPair(mMeanLocation, mPrecisionShape / mPrecisionInvScale);
                }
            }
        }

        /// <summary>
        /// Evaluates the probability density function for a NormalGamma distribution.
        /// </summary>
        public double Density(MeanPrecisionPair mp)
        {
            return Density(mp.Mean, mp.Precision);
        }

        /// <summary>
        /// Evaluates the probability density function for a NormalGamma distribution.
        /// </summary>
        public double Density(double mean, double prec)
        {
            if (Double.IsPositiveInfinity(mPrecisionInvScale) && mMeanScale == 0.0)
            {
                throw new NotImplementedException();
            }
            else if (Double.IsPositiveInfinity(mPrecisionInvScale))
            {
                throw new NotImplementedException();
            }
            else if (mMeanScale == 0.0)
            {
                throw new NotImplementedException();
            }
            else
            {
                double e = -0.5 * prec * (mean - mMeanLocation) * (mean - mMeanLocation) - prec * mPrecisionInvScale;
                return System.Math.Pow(prec * mPrecisionInvScale, mPrecisionShape) * System.Math.Exp(e)
                        / (Math.Constants.Sqrt2Pi * System.Math.Sqrt(prec) * Math.SpecialFunctions.Gamma(mPrecisionShape));
            }
        }

        /// <summary>
        /// Evaluates the log probability density function for a NormalGamma distribution.
        /// </summary>
        public double DensityLn(MeanPrecisionPair mp)
        {
            return DensityLn(mp.Mean, mp.Precision);
        }

        /// <summary>
        /// Evaluates the log probability density function for a NormalGamma distribution.
        /// </summary>
        public double DensityLn(double mean, double prec)
        {
            if (Double.IsPositiveInfinity(mPrecisionInvScale) && mMeanScale == 0.0)
            {
                throw new NotImplementedException();
            }
            else if (Double.IsPositiveInfinity(mPrecisionInvScale))
            {
                throw new NotImplementedException();
            }
            else if (mMeanScale == 0.0)
            {
                throw new NotImplementedException();
            }
            else
            {
                double e = -0.5 * prec * (mean - mMeanLocation) * (mean - mMeanLocation) - prec * mPrecisionInvScale;
                return (mPrecisionShape - 0.5) * System.Math.Log(prec) + mPrecisionShape * System.Math.Log(mPrecisionInvScale) + e
                        - Math.Constants.LogSqrt2Pi - Math.SpecialFunctions.GammaLn(mPrecisionShape);
            }
        }

        /// <summary>
        /// Samples a NormalGamma distributed random variable.
        /// </summary>
        /// <returns>A random number from this distribution.</returns>
        public MeanPrecisionPair Sample()
        {
            return NormalGamma.Sample(RandomNumberGenerator, mMeanLocation, mMeanScale, mPrecisionShape, mPrecisionInvScale);
        }

        /// <summary>
        /// Samples an array of NormalGamma distributed random variables.
        /// </summary>
        /// <param name="size">The number of variables needed.</param>
        /// <returns>An array of random numbers from this distribution.</returns>
        public MeanPrecisionPair[] Sample(int size)
        {
            return NormalGamma.Sample(RandomNumberGenerator, size, mMeanLocation, mMeanScale, mPrecisionShape, mPrecisionInvScale);
        }

        /// <summary>
        /// Checks the parameters of a NormalGamma distribution.
        /// </summary>
        /// <param name="meanScale">The scale of the mean.</param>
        /// <param name="precShape">The shape of the precision.</param>
        /// <param name="precInvScale">The inverse scale of the precision.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the mean scale is negative.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the inverse precision scale is negative.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the precision shape is negative.</exception>
        private static void CheckParameters(double meanScale, double precShape, double precInvScale)
        {
            if (meanScale < 0.0)
            {
                throw new ArgumentOutOfRangeException("meanScale", Resources.ParameterCannotBeNegative);
            }
            else if (precShape <= 0.0)
            {
                throw new ArgumentOutOfRangeException("precShape", Resources.ParameterCannotBeNegative);
            }
            else if (precInvScale <= 0.0)
            {
                throw new ArgumentOutOfRangeException("precInvScale", Resources.ParameterCannotBeNegative);
            }
        }

        /// <summary>
        /// Samples an array of NormalGamma distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="meanLocation">The location of the mean.</param>
        /// <param name="meanScale">The scale of the mean.</param>
        /// <param name="precShape">The shape of the precision.</param>
        /// <param name="precInvScale">The inverse scale of the precision.</param>
        public static MeanPrecisionPair Sample(System.Random rnd, double meanLocation, double meanScale, double precShape, double precInvScale)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(meanScale, precShape, precInvScale);
            }

            MeanPrecisionPair mp = new MeanPrecisionPair();

            // Sample the precision.
            if(Double.IsPositiveInfinity(precInvScale))
            {
                mp.Precision = precShape;
            }
            else
            {
                mp.Precision = Gamma.Sample(rnd, precShape, precInvScale);
            }

            // Sample the mean.
            if (meanScale == 0.0)
            {
                mp.Mean = meanLocation;
            }
            else
            {
                mp.Mean = Normal.Sample(rnd, meanLocation, System.Math.Sqrt(meanScale / mp.Precision));
            }

            return mp;
        }

        /// <summary>
        /// Samples an array of NormalGamma distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="meanLocation">The location of the mean.</param>
        /// <param name="meanScale">The scale of the mean.</param>
        /// <param name="precShape">The shape of the precision.</param>
        /// <param name="precInvScale">The inverse scale of the precision.</param>
        public static MeanPrecisionPair[] Sample(System.Random rnd, int n, double meanLocation, double meanScale, double precShape, double precInvScale)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(meanScale, precShape, precInvScale);
            }

            // First sample all the precisions independently.
            double[] precs = null;
            if (Double.IsPositiveInfinity(precInvScale))
            {
                precs = new double[n];
                for (int i = 0; i < n; i++)
                {
                    precs[i] = precShape;
                }
            }
            else
            {
                precs = Gamma.Sample(rnd, n, precShape, precInvScale);
            }

            // Construct all the mean precision pairs.
            MeanPrecisionPair[] arr = new MeanPrecisionPair[n];

            // Conditionally sample all the mean.
            for (int i = 0; i < n; i++)
            {
                arr[i].Precision = precs[i];
                if (meanScale == 0.0)
                {
                    arr[i].Mean = meanLocation;
                }
                else
                {
                    arr[i].Mean = Normal.Sample(rnd, meanLocation, System.Math.Sqrt(meanScale / precs[i]));
                }
            }

            return arr;
        }
    }
}